$title	ORANI-G: A Template Implementation with MPSGE (Top-Down)

$if not set dataset $set dataset ozdat934

$call 'har2gdx %dataset%'

set	c(*)	Commodities,
	i(*)	Industries,
	l(*)	Labor occupations,
	r(*)	Regions
	s(*)	Sources ;

$gdxin %dataset%
$load c=com i=ind l=occ r=reg s=src


set	mrg(c)	Margin commodities;

$load mrg=mar

PARAMETER
	v1bas(c,s,i)		"Intermediate Basic [[1BAS:COM*SRC*IND]]",
	v2bas(c,s,i)		"Investment Basic [[2BAS:COM*SRC*IND]]"
	v3bas(c,s)		"Households Basic [[3BAS:COM*SRC]]"
	v4bas(c)		"Exports Basic [[4BAS:COM]]"
	v5bas(c,s)		"Government Basic [[5BAS:COM*SRC]]"
	v6bas(c,s)		"Inventory Changes [[6BAS:COM*SRC]]"
	v1mar(c,s,i,mrg)	"Intermediate Margins [[1MAR:COM*SRC*IND*MAR]]"
	v2mar(c,s,i,mrg)	"Investment Margins [[2MAR:COM*SRC*IND*MAR]]"
	v3mar(c,s,mrg)		"Households Margins [[3MAR:COM*SRC*MAR]]"
	v4mar(c,mrg)		"Exports Margins [[4MAR:COM*MAR]]"
	v5mar(c,s,mrg)		"Government Margins [[5MAR:COM*SRC*MAR]]  "
	v1tax(c,s,i)		"Intermediate Tax [[1TAX:COM*SRC*IND]]"
	v2tax(c,s,i)		"Investment Tax [[2TAX:COM*SRC*IND]]"
	v3tax(c,s)		"Households Tax [[3TAX:COM*SRC]]"
	v4tax(c)		"Exports Tax [[4TAX:COM]]"
	v5tax(c,s)		"Government Tax [[5TAX:COM*SRC]]"
	v1cap(i)		"Capital [[1CAP:IND]]"
	v1lab(i,l)		"Labour [[1LAB:IND*OCC]]"
	v1lnd(i)		"Land [[1LND:IND]]"
	v1ptx(i)		"Production Tax [[1PTX:IND]]  "
	v1oct(i)		"Other Costs [[1OCT:IND]]"
	make(c,i)		"Multiproduct Matrix [[MAKE:COM*IND]]"
	v0tar(c)		"Tariff Revenue [[0TAR:COM]]"

	sigma1lab(i)		"Labour Sigma [[SLAB:IND]]"
	sigma1prim(i)		"Primary Factor Sigma [[P028:IND]]"
	sigma1(c)		"Intermediate Armington [[1ARM:COM]]"
	sigma1out(i)		"Output Sigma [[SCET:IND]]"
	sigma2(c)		"Investment Armington [[2ARM:COM]]"
	sigma3(c)		"Final demand Argminton",
	exp_elast(c)		"Individual Export Elasticities [[P018:COM]]"
	isindivexp(*)		"'Flag, >0.5 for individual export commodities, else collective export [[ITEX:COM]]'"
	isloccom(c)		"'Flag for regional extension, >0.5 for local coms, else national [[LCOM:COM]]'"

	regshr1(i,r)		"Regional output shares [[R001:IND*REG]]"
	regshr2(i,r)		"Regional investment shares [[R002:IND*REG]]"
	regshr3(c,r)		"Regional consumption shares [[R003:COM*REG]]"
	regshr4(c,r)		"Regional export shares [[R004:COM*REG]]"
	regshr5(c,r)		"Regional government shares [[R005:COM*REG]]"
	regshr6(c,r)		"Regional inventory shares [[R006:COM*REG]]";

scalar	frisch			"Frisch  [[P021:]]";

$load v1bas v2bas v3bas v4bas v5bas v6bas v1mar v2mar v3mar v4mar v5mar v1tax v2tax v3tax 
$load v4tax v5tax v1cap v1lab v1lnd v1ptx v1oct make v0tar sigma1lab sigma1prim sigma1 
$load sigma1out sigma2 sigma3 frisch exp_elast isindivexp isloccom regshr1 
$load regshr2 regshr3 regshr4 regshr5 regshr6 
$load

parameter	regchk		Check sums on regional totals;

regchk(i,"1") = sum(r, regshr1(i,r));
regchk(i,"2") = sum(r, regshr2(i,r));
regchk(c,"3") = sum(r, regshr3(c,r));
regchk(c,"4") = sum(r, regshr4(c,r));
regchk(c,"5") = sum(r, regshr5(c,r));
regchk(c,"6") = sum(r, regshr6(c,r));

*	Filter the data to avoid spurious entries:

v1bas(c,s,i)$(not round(v1bas(c,s,i), 6)) = 0;
v2bas(c,s,i)$(not round(v2bas(c,s,i), 6)) = 0;
v3bas(c,s)$(not round(v3bas(c,s), 6)) = 0;
v4bas(c)$(not round(v4bas(c), 6)) = 0;
v5bas(c,s)$(not round(v5bas(c,s), 6)) = 0;
v6bas(c,s)$(not round(v6bas(c,s), 6)) = 0;
v1mar(c,s,i,mrg)$(not round(v1mar(c,s,i,mrg), 6)) = 0;
v2mar(c,s,i,mrg)$(not round(v2mar(c,s,i,mrg), 6)) = 0;
v3mar(c,s,mrg)$(not round(v3mar(c,s,mrg), 6)) = 0;
v4mar(c,mrg)$(not round(v4mar(c,mrg), 6)) = 0;
v5mar(c,s,mrg)$(not round(v5mar(c,s,mrg), 6)) = 0;
v1tax(c,s,i)$(not round(v1tax(c,s,i), 6)) = 0;
v2tax(c,s,i)$(not round(v2tax(c,s,i), 6)) = 0;
v3tax(c,s)$(not round(v3tax(c,s), 6)) = 0;
v4tax(c)$(not round(v4tax(c), 6)) = 0;
v5tax(c,s)$(not round(v5tax(c,s), 6)) = 0;
v1cap(i)$(not round(v1cap(i), 6)) = 0;
v1lab(i,l)$(not round(v1lab(i,l), 6)) = 0;
v1lnd(i)$(not round(v1lnd(i), 6)) = 0;
v1ptx(i)$(not round(v1ptx(i), 6)) = 0;
v1oct(i)$(not round(v1oct(i), 6)) = 0;
make(c,i)$(not round(make(c,i), 6)) = 0;
v0tar(c)$(not round(v0tar(c), 6)) = 0;
sigma1lab(i)$(not round(sigma1lab(i), 6)) = 0;
sigma1prim(i)$(not round(sigma1prim(i), 6)) = 0;
sigma1(c)$(not round(sigma1(c), 6)) = 0;
sigma1out(i)$(not round(sigma1out(i), 6)) = 0;
sigma2(c)$(not round(sigma2(c), 6)) = 0;


parameter	scalefac	Scale factor /1000/;

v1bas(c,s,i) = v1bas(c,s,i)/scalefac;
v2bas(c,s,i) = v2bas(c,s,i)/scalefac;
v3bas(c,s) = v3bas(c,s)/scalefac;
v4bas(c) = v4bas(c)/scalefac;
v5bas(c,s) = v5bas(c,s)/scalefac;
v6bas(c,s) = v6bas(c,s)/scalefac;
v1mar(c,s,i,mrg) = v1mar(c,s,i,mrg)/scalefac;
v2mar(c,s,i,mrg) = v2mar(c,s,i,mrg)/scalefac;
v3mar(c,s,mrg) = v3mar(c,s,mrg)/scalefac;
v4mar(c,mrg) = v4mar(c,mrg)/scalefac;
v5mar(c,s,mrg) = v5mar(c,s,mrg)/scalefac;
v1tax(c,s,i) = v1tax(c,s,i)/scalefac;
v2tax(c,s,i) = v2tax(c,s,i)/scalefac;
v3tax(c,s) = v3tax(c,s)/scalefac;
v4tax(c) = v4tax(c)/scalefac;
v5tax(c,s) = v5tax(c,s)/scalefac;
v1cap(i) = v1cap(i)/scalefac;
v1lab(i,l) = v1lab(i,l)/scalefac;
v1lnd(i) = v1lnd(i)/scalefac;
v1ptx(i) = v1ptx(i)/scalefac;
v1oct(i) = v1oct(i)/scalefac;
make(c,i) = make(c,i)/scalefac;
v0tar(c) = v0tar(c)/scalefac;


parameter
	basicflow(c,s)		Basic flows
	marginflow(c,s,mrg)	Margin flows;

basicflow(c,s) = sum(i, v1bas(c,s,i) + v2bas(c,s,i)) + v3bas(c,s) + v4bas(c)$sameas(s,"dom") 
				+ v5bas(c,s) + v6bas(c,s);
marginflow(c,s,mrg) = sum(i, v1mar(c,s,i,mrg) + v2mar(c,s,i,mrg)) + v3mar(c,s,mrg) 
				+ v4mar(c,mrg)$sameas(s,"dom") + v5mar(c,s,mrg);

parameter	comchk(c,*)	Commodity flow check;
comchk(c,"supply") = sum(i,make(c,i));
comchk(c,"basic")=  basicflow(c,"dom");
comchk(mrg,"margin")= sum((c,s), marginflow(c,s,mrg));
comchk(c,"chksum") = comchk(c,"supply") - comchk(c,"basic") - comchk(c,"margin");
display comchk;

set cst/use,margin,labor,capital,land,taxes,prodtax,othercost, profit/;
parameter	vshare(i,cst)	Benchmark value shares;
vshare(i,"use") = sum((c,s), v1bas(c,s,i));
vshare(i,"margin") = sum((c,s,mrg), v1mar(c,s,i,mrg));
vshare(i,"labor") = sum(l,v1lab(i,l));
vshare(i,"capital") = v1cap(i); 
vshare(i,"land") = v1lnd(i);
vshare(i,"taxes") = sum((c,s), v1tax(c,s,i));
vshare(i,"prodtax") = v1ptx(i);
vshare(i,"othercost") = v1oct(i);
vshare(i,"profit") = sum(cst,vshare(i,cst)) - sum(c, make(c,i));
display vshare;

parameter	t1(c,s,i)	Tax rates on intermediate demand ,
		t2(c,s,i)	Tax rates on investment demand,
		t3(c,s)		Tax rates on final demand,
		t4(c)		Tax rates on exports, 
		t5(c,s)		Tax rates on government demand, 
		v1(c,i)		Armington composite intermediate demand,
		v2(c,i)		Armington composite investment demand,
		v3(c)		Armington composite household demand,
		v4(c)		Exports,
		v5(c)		Armington government demand, 
		vinv(i)		Aggregate investment,
		vgovt		Aggregate public demand, 
		vbopdef		Base year current account, 
		vm(c)		Imports,
		vlstax		Implicit direct taxes, 
		vmake(c)	Aggregate output,
		vlab(l)		Value of labor income,
		vlnd		Value of land ents,
		vcap		Value of capital rent, 
		vhhd(r)		Value of household demand;

t1(c,s,i)$v1tax(c,s,i) = v1tax(c,s,i) / (v1bas(c,s,i)+sum(mrg, v1mar(c,s,i,mrg)));
vlab(l) = sum(i,v1lab(i,l));
v1cap(i) = v1cap(i)+min(0,v1oct(i));
v1oct(i) = max(0, v1oct(i));
vcap = sum(i, v1cap(i));
vlnd = sum(i, v1lnd(i));
vmake(c) = sum(i, make(c,i));

v1(c,i) = sum(s, v1bas(c,s,i) + v1tax(c,s,i) + sum(mrg, v1mar(c,s,i,mrg) ));
v2(c,i) = sum(s, v2bas(c,s,i) + v2tax(c,s,i) + sum(mrg, v2mar(c,s,i,mrg)));
t2(c,s,i)$v2tax(c,s,i) = v2tax(c,s,i) /(v2bas(c,s,i) + sum(mrg, v2mar(c,s,i,mrg)));
v4(c) = sum(mrg, v4mar(c,mrg)) + v4bas(c) + v4tax(c);
v3(c)= sum(s, v3bas(c,s) + v3tax(c,s) + sum(mrg, v3mar(c,s,mrg)));
vhhd(r) = sum(c, v3(c)*regshr3(c,r));

t3(c,s)$v3tax(c,s) = v3tax(c,s) / (v3bas(c,s) + sum(mrg, v3mar(c,s,mrg)));
vinv(i) = sum(c, v2(c,i));
t4(c)$v4(c) = v4tax(c)/v4(c);
v5bas(c,s) = max(0,v5bas(c,s));
v5(c) = sum(s, v5bas(c,s) + v5tax(c,s) + sum(mrg, v5mar(c,s,mrg)));
t5(c,s)$v5tax(c,s) = v5tax(c,s) / (v5bas(c,s) + sum(mrg, v5mar(c,s,mrg)));
vgovt = sum(c, v5(c));
vm(c) = sum(i, v1bas(c,"imp",i)+v2bas(c,"imp",i)) + v3bas(c,"imp") + v5bas(c,"imp");
vbopdef = sum(c, vm(c) - v4(c));
vlstax = vgovt - (sum((c,s,i), v1tax(c,s,i)+v2tax(c,s,i)) + sum((c,s), v3tax(c,s))
	+ sum(c, v4tax(c)) + sum((c,s),v5tax(c,s)) + vbopdef);

parameter vrtax;
vrtax("total") =   sum(c,    v3(c)    )
	   + sum(c,    v6bas(c,"dom") )
	   + sum(i,    vinv(i)        )
	   - sum((i,l),v1lab(i,l)     )
	   - sum(i,    v1lnd(i)       )
	   - sum(i,    v1cap(i)       )
	   - sum(i,    v1oct(i)       );

vrtax(r) =   sum(c,    v3(c)          * regshr3(c,r))
	   + sum(c,    v6bas(c,"dom") * regshr6(c,r)) 
	   + sum(i,    vinv(i)        * regshr2(i,r))
	   - sum((i,l),v1lab(i,l)     * regshr1(i,r))
	   - sum(i,    v1lnd(i)       * regshr1(i,r))
	   - sum(i,    v1cap(i)       * regshr1(i,r))
	   - sum(i,    v1oct(i)       * regshr1(i,r));

vrtax("totchk") = vrtax("total")   + vlstax;
vrtax("regchk") = sum(r, vrtax(r)) + vlstax;
display vrtax;


parameter p1ref, p2ref, p3ref, p5ref;
p1ref(c,s,i) = 1 + t1(c,s,i);
p2ref(c,s,i) = 1 + t2(c,s,i);
p3ref(c,s) = 1 + t3(c,s);
p5ref(c,s) = 1 + t5(c,s);

parameter dmkt;
dmkt(c,"supply") = vmake(c);
dmkt(c,"basic")   = sum(i, v1bas(c,"dom",i)+v2bas(c,"dom",i)) + v4bas(c) 
			+ v3bas(c,"dom") + v5bas(c,"dom") + v6bas(c,"dom");
dmkt(mrg,"margin") = sum((c,s,i), v1mar(c,s,i,mrg)+v2mar(c,s,i,mrg))
			+ sum((c,s),v3mar(c,s,mrg)) + sum(c, v4mar(c,mrg))
			+ sum((c,s),v5mar(c,s,mrg));
dmkt(c,"chk")= dmkt(c,"supply") - dmkt(c,"basic") - dmkt(c,"margin");
display dmkt;

$ontext

$model:oranig

$sectors:

	hhd		! Household demand
	govd		! Government demand

	y(i)		! Production
	invd(i)		! Investment demand

	e(c)$v4(c)	! Exports
	m(c)$vm(c)	! Imports

	x1(c,i)$v1(c,i)	! Armington composite inputs to production
	x2(c,i)$v2(c,i)	! Armington composite inputs to investment
	x3(c)$v3(c)	! Armington composite inputs to demand
	x4(c)		! Armington composite of exports and domestic
	x5(c)		! Armington composite inputs to government

$commodities:

	phh		! Household demand price 
	pgov		! Government demand price 
	pinv(i)$vinv(i)	! Investment demand price 
	pfx		! Foreign exchange

	pd(c)		! Domestic price
	pm(c)$vm(c)	! Import price

	pn		! Land price
	rk		! Capital return
	py(c)		! Output prices
	pl(l)		! Wage rate
	ps(i)$v1oct(i)	! Other cost (specific input)

	p1(c,i)$v1(c,i)	! Armington composite (production)
	p2(c,i)$v2(c,i)	! Armington composite (investment)
	p3(c)$v3(c)	! Armington composite (final demand)
	p4(c)$v4bas(c)	! Domestic supply price
	p5(c)$v5(c)	! Armington composite (government)

$consumers:
	hh		! Household
	govt		! Government

$auxiliary:
	lstax		! Lumpsum tax replacement

*	-------------------------------------------------------------------

$demand:hh
	e:pd(c)		q:(-v6bas(c,"dom"))
	e:pinv(i)	q:(-vinv(i))
	e:pl(l)		q:vlab(l)
	e:pn		q:vlnd
	e:rk		q:vcap
	e:ps(i)		q:v1oct(i)
	e:pgov		q:(-vlstax)
	e:pfx		q:(-sum(r,vhhd(r)))		r:lstax
	d:phh

$demand:govt
	e:pfx		q:(sum(r,vhhd(r)))	r:lstax
	e:pfx		q:vbopdef
	e:pgov		q:vlstax
	d:pgov		

*	-------------------------------------------------------------------

$prod:m(c)$vm(c)  s:0
	o:pm(c)		q:vm(c)
	i:pfx		q:vm(c)

*	-------------------------------------------------------------------

$prod:y(i)	t:sigma1out(i)  s:0  m:0 va:sigma1prim(i) ld(va):sigma1lab(i)
	o:py(c)		q:make(c,i)
	i:pn		q:v1lnd(i)	va:
	i:pl(l)		q:v1lab(i,l)	ld:
	i:rk		q:v1cap(i)	va:
	i:ps(i)		q:v1oct(i)	va:
	i:p1(c,i)	q:v1(c,i)	m:

$prod:x1(c,i)$v1(c,i)  s:sigma1(c)  d:0  m:0
	o:p1(c,i)	q:v1(c,i)
	i:pd(c)		q:v1bas(c,"dom",i)	a:govt	t:t1(c,"dom",i)  d: p:p1ref(c,"dom",i)
	i:pd(mrg)	q:v1mar(c,"dom",i,mrg)		t:t1(c,"dom",i)  d: p:p1ref(c,"dom",i)
	i:pm(c)		q:v1bas(c,"imp",i)		t:t1(c,"imp",i)  m: p:p1ref(c,"imp",i)
	i:pd(mrg)	q:v1mar(c,"imp",i,mrg)		t:t1(c,"imp",i)  m: p:p1ref(c,"imp",i)


*	-------------------------------------------------------------------

$prod:invd(i) s:1
	o:pinv(i)	q:vinv(i)
	i:p2(c,i)	q:v2(c,i)

$prod:x2(c,i)$v2(c,i)  s:sigma2(c) d:0  m:0
	o:p2(c,i)	q:v2(c,i)
	i:pd(c)		q:v2bas(c,"dom",i)	a:govt	t:t2(c,"dom",i)  d: p:p2ref(c,"dom",i)
	i:pd(mrg)	q:v2mar(c,"dom",i,mrg)		t:t2(c,"dom",i)  d: p:p2ref(c,"dom",i)
	i:pm(c)		q:v2bas(c,"imp",i)		t:t2(c,"imp",i)  m: p:p2ref(c,"imp",i)
	i:pd(mrg)	q:v2mar(c,"imp",i,mrg)		t:t2(c,"imp",i)  m: p:p2ref(c,"imp",i)


*	-------------------------------------------------------------------

$prod:hhd s:1
	o:phh		q:(sum(r,vhhd(r)))
	i:p3(c)		q:v3(c)

$prod:x3(c)$v3(c) s:sigma3(c) d:0  m:0
	o:p3(c)		q:v3(c)
	i:pd(c)		q:v3bas(c,"dom")	a:govt	t:t3(c,"dom")  d: p:p3ref(c,"dom")
	i:pd(mrg)	q:v3mar(c,"dom",mrg)		t:t3(c,"dom")  d: p:p3ref(c,"dom")
	i:pm(c)		q:v3bas(c,"imp")		t:t3(c,"imp")  m: p:p3ref(c,"imp")
	i:pd(mrg)	q:v3mar(c,"imp",mrg)		t:t3(c,"imp")  m: p:p3ref(c,"imp")

*	-------------------------------------------------------------------

$prod:e(c)$v4(c)   s:0
	o:pfx		q:v4(c)			a:govt	t:t4(c)
	i:pd(mrg)	q:v4mar(c,mrg)
	i:p4(c)		q:v4bas(c)

$prod:x4(c)   t:2
	o:pd(c)		q:(vmake(c)-v4bas(c))
	o:p4(c)		q:v4bas(c)
	i:py(c)		q:vmake(c)

*	-------------------------------------------------------------------

$prod:govd
	o:pgov		q:vgovt
	i:p5(c)		q:v5(c)

$prod:x5(c)$v5(c) s:sigma3(c) d:0  m:0
	o:p5(c)		q:v5(c)
	i:pd(c)		q:v5bas(c,"dom")	a:govt	t:t5(c,"dom")  d: p:p5ref(c,"dom")
	i:pd(mrg)	q:v5mar(c,"dom",mrg)		t:t5(c,"dom")  d: p:p5ref(c,"dom")
	i:pm(c)		q:v5bas(c,"imp")		t:t5(c,"imp")  m: p:p5ref(c,"imp")
	i:pd(mrg)	q:v5mar(c,"imp",mrg)		t:t5(c,"imp")  m: p:p5ref(c,"imp")

$constraint:lstax
	govd =e= 1;

$offtext
$sysinclude mpsgeset oranig

oranig.iterlim = 0;
$include oranig.gen
solve oranig using mcp;

parameter	income(*)	Income level by region;

income(r) = 
	sum(c, pd.l(c)		* (-v6bas(c,"dom")*regshr6(c,r)))
	+ sum(i, pinv.l(i)	* (-vinv(i)*regshr2(i,r)))
	+ sum(l, pl.l(l)	* (sum(i,v1lab(i,l)*regshr1(i,r))))
	+ pn.l			* (sum(i, v1lnd(i)*regshr1(i,r)))
	+ rk.l			* (sum(i, v1cap(i)*regshr1(i,r)))
	+ sum(i, ps.l(i)	* (v1oct(i)*regshr1(i,r)))
	+ pgov.l		* vrtax(r)
	+ pfx.l		* (-vhhd(r))	* lstax.l;


income("chk") = sum(r, income(r)) - hh.l;

parameter	ev(*,*)	Regional equilvalent variation in income;
ev(r,"%") = 100 * (income(r)/(phh.l*vhhd(r))-1);
display income, ev;

*	Consider a free-trade calculation with lumpsum replacement:

t1(c,"imp",i) = 0;
t2(c,"imp",i) = 0;
t3(c,"imp") = 0;
t5(c,"imp") = 0;

oranig.iterlim = 10000;
$include oranig.gen
solve oranig using mcp;

income(r) = 
	sum(c, pd.l(c)		* (-v6bas(c,"dom")*regshr6(c,r)))
	+ sum(i, pinv.l(i)	* (-vinv(i)*regshr2(i,r)))
	+ sum(l, pl.l(l)	* (sum(i,v1lab(i,l)*regshr1(i,r))))
	+ pn.l			* (sum(i, v1lnd(i)*regshr1(i,r)))
	+ rk.l			* (sum(i, v1cap(i)*regshr1(i,r)))
	+ sum(i, ps.l(i)	* (v1oct(i)*regshr1(i,r)))
	+ pgov.l		* vrtax(r)
	+ pfx.l		* (-vhhd(r))	* lstax.l;


income("chk") = sum(r, income(r)) - hh.l;
display income;
ev(r,"%") = 100 * (income(r)/(phh.l*vhhd(r))-1);
ev(r,"value") = income(r)/phh.l-vhhd(r);
ev("total","value") =sum(r, income(r)/phh.l-vhhd(r));
display ev;
